PCA 可视化 与 DimReduc 类




() DimReduc 类的定义
./seurat-object-4.0.4/R/dimreduc.R:36:DimReduc <- setClass(

每个降维类包含：
- a cell embeddings matrix, (必须)
- a feature loadings matrix, (可选)
- and a projected feature loadings matrix. (可选)


#' The Dimensional Reduction Class
#'
#' The DimReduc object stores a dimensionality reduction taken out in Seurat;
#' each DimReduc consists of a cell embeddings matrix, a feature loadings
#' matrix, and a projected feature loadings matrix.
#'
#' @slot cell.embeddings Cell embeddings matrix (required)
#' @slot feature.loadings Feature loadings matrix (optional)
#' @slot feature.loadings.projected Projected feature loadings matrix (optional)
#' @slot assay.used Name of assay used to generate \code{DimReduc} object
#' @slot global Is this \code{DimReduc} global/persistent? If so, it will not be
#' removed when removing its associated assay
#' @slot stdev A vector of standard deviations
#' @slot key Key for the \code{DimReduc}, must be alphanumeric characters
#' followed by an underscore
#' @slot jackstraw A \code{\link{JackStrawData-class}} object associated with
#' this \code{DimReduc}
#' @slot misc Utility slot for storing additional data associated with the
#' \code{DimReduc} (e.g. the total variance of the PCA)
#'
#' @name DimReduc-class
#' @rdname DimReduc-class
#' @exportClass DimReduc
#'
DimReduc <- setClass(
  Class = 'DimReduc',
  slots = c( #每个slot的用途，我的猜测，不一定准确
    cell.embeddings = 'matrix', #细胞的位置
    feature.loadings = 'matrix', #转置矩阵
    feature.loadings.projected = 'matrix', #不确定
    assay.used = 'character', #实验名字, assays list中的name，"RNA"
    global = 'logical', #如果是全局的，则不能轻易删除
    stdev = 'numeric', #PCA的碎石图就是用的这个做y轴
    key = 'character', #前缀 "PC_"
    jackstraw = 'JackStrawData', #模拟数据类，又一层内部类
    misc = 'list' #杂项，是一个list()
  )
)











() CreateDimReducObject()
./seurat-object-4.0.4/R/dimreduc.R:88:CreateDimReducObject <- function(

#' Create a DimReduc object
#'
#' @param embeddings A matrix with the cell embeddings
#' @param loadings A matrix with the feature loadings
#' @param projected A matrix with the projected feature loadings
#' @param assay Assay used to calculate this dimensional reduction
#' @param stdev Standard deviation (if applicable) for the dimensional reduction
#' @param key A character string to facilitate looking up features from a
#' specific DimReduc
#' @param global Specify this as a global reduction (useful for visualizations)
#' @param jackstraw Results from the JackStraw function
#' @param misc list for the user to store any additional information associated
#' with the dimensional reduction
#'
#' @return A \code{\link{DimReduc}} object
#'
#' @aliases SetDimReduction
#'
#' @export
#'
#' @concept dimreduc
#'
#' @examples
#' data <- GetAssayData(pbmc_small[["RNA"]], slot = "scale.data")
#' pcs <- prcomp(x = data)
#' pca.dr <- CreateDimReducObject(
#'   embeddings = pcs$rotation,
#'   loadings = pcs$x,
#'   stdev = pcs$sdev,
#'   key = "PC",
#'   assay = "RNA"
#' )
#'
CreateDimReducObject <- function(
  embeddings = new(Class = 'matrix'),
  loadings = new(Class = 'matrix'),
  projected = new(Class = 'matrix'),
  assay = NULL,
  stdev = numeric(),
  key = NULL,
  global = FALSE,
  jackstraw = NULL,
  misc = list()
) {
  # (A1) 如果 assay 为空，则给出警告，并默认用 "RNA"，就是 pbmc_small@assays$RNA
  if (is.null(x = assay)) {
    warning(
      "No assay specified, setting assay as RNA by default.",
      call. = FALSE,
      immediate. = TRUE
    )
    assay <- "RNA"
  }

  # (A2) 如果 key 为空，且 embeddings 列名为空，则报错
  # Try to infer key from column names
  if (is.null(x = key) && is.null(x = colnames(x = embeddings))) {
    stop("Please specify a key for the DimReduc object")
  # (B2)如果 key 为空，则自动匹配 字母数字 到下划线部分，并取唯一值
  } else if (is.null(x = key)) {
    key <- regmatches(
      x = colnames(x = embeddings),
      m = regexec(pattern = '^[[:alnum:]]+_', text = colnames(x = embeddings))
    )
    key <- unique(x = unlist(x = key, use.names = FALSE)) #"PC_"
  }

  # (A3)如果key长度不等于1，则报错(主要是针对 A2-B2部分做验证)
  if (length(x = key) != 1) {
    stop("Please specify a key for the DimReduc object")
  # 如果key长度为1，且不符合格式: 字母数字若干 下划线结尾
  } else if (!grepl(pattern = '^[[:alnum:]]+_$', x = key)) {
    old.key  <- key
    key <- UpdateKey(key = old.key) #重新更新 key(匹配到 字母数字的，后面加下划线；只有下划线的，生成3个随机字母)，见下文
    # 替换列名中的key子字符串
    colnames(x = embeddings) <- gsub(
      x = colnames(x = embeddings),
      pattern = old.key,
      replacement = key
    )
    #因为程序纠正过格式，所以给出警告：key必须 一个或多个字母数字，然后下划线结尾
    warning(
      "All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to ",
      key,
      call. = FALSE,
      immediate. = TRUE
    )
  }

  # (A4) 保证 embeddings 的后缀是数字
  # ensure colnames of the embeddings are the key followed by a numeric
  # 如果没有列名，则警告并自动添加：key 下划线 递增的数字
  if (is.null(x = colnames(x = embeddings))) {
    warning(
      "No columnames present in cell embeddings, setting to '",
      key,
      "1:",
      ncol(x = embeddings),
      "'",
      call. = FALSE,
      immediate. = TRUE
    )
    colnames(x = embeddings) <- paste0(key, 1:ncol(x = embeddings))
  # 如果有列名，且不是都符合格式： key开头，后面是数字结尾
  } else if (!all(grepl(pattern = paste0('^', key, "[[:digit:]]+$"), x = colnames(x = embeddings)))) {
    # 退一步，匹配数字结尾，比如 x=c("PC_3", "tSNE_2")，返回的是 "3" "2"
    digits <- unlist(x = regmatches(
      x = colnames(x = embeddings),
      m = regexec(pattern = '[[:digit:]]+$', text = colnames(x = embeddings))
    ))
    # 如果长度不等，则报错。也即是有些不是数字结尾。
    if (length(x = digits) != ncol(x = embeddings)) {
      stop("Please ensure all column names in the embeddings matrix are the key plus a digit representing a dimension number")
    }
    # 拼接为 key 后面加数字的形式
    colnames(x = embeddings) <- paste0(key, digits)
  }


  # (A5) 如果 loadings 不是空矩阵
  if (!IsMatrixEmpty(x = loadings)) { #源码解析 13-2.9 讲过 IsMatrixEmpty()
    # 如果行名有空的，则报错
    if (any(rownames(x = loadings) == '')) {
      stop("Feature names of loadings matrix cannot be empty", call. = FALSE)
    }
    # gene loadings 列名 等于 cell embedding 列名
    colnames(x = loadings) <- colnames(x = embeddings)
  }

  # (A6) 如果 projected 不是空矩阵
  if (!IsMatrixEmpty(x = projected)) {
    # 如果行名有空的，则报错
    if (any(rownames(x = loadings) == '')) {
      stop("Feature names of projected loadings matrix cannot be empty", call. = FALSE)
    }
    # 投影列名 =  embedding 列名。
    colnames(x = projected) <- colnames(x = embeddings)
  }

  #(A7) jackstraw 如果为空，则新建一个空对象。
  jackstraw <- jackstraw %||% new(Class = 'JackStrawData')
  
  #(A8) 创建 降维 对象
  dim.reduc <- new(
    Class = 'DimReduc',
    cell.embeddings = embeddings,
    feature.loadings = loadings,
    feature.loadings.projected = projected,
    assay.used = assay,
    global = global,
    stdev = stdev,
    key = key,
    jackstraw = jackstraw,
    misc = misc
  )
  return(dim.reduc)
}










() print.DimReduc()
./seurat-object-4.0.4/R/dimreduc.R:611:print.DimReduc <- function(


#' @describeIn DimReduc-methods Prints a set of features that most strongly
#' define a set of components; \strong{note}: requires feature loadings to be
#' present in order to work  #依赖 feature loadings
#'
#' @param dims Number of dimensions to display
#' @param nfeatures Number of genes to display
#' @param projected Use projected slot
#' @param ... Arguments passed to other methods
#'
#' @return \code{print}: Displays set of features defining the components and
#' invisibly returns \code{x}
#'
#' @aliases print
#' @seealso \code{\link[base]{cat}}
#'
#' @export
#' @method print DimReduc
#'
print.DimReduc <- function(
  x,  #第一个参数，有时候是 x，有时候是 object，有规律吗？
  dims = 1:5,
  nfeatures = 20,
  projected = FALSE,
  ...
) {
  CheckDots(...)
  # (A1) 就是 pbmc_small@reductions$pca@feature.loadings，每个基因的权重
  loadings <- Loadings(object = x, projected = projected) #源码解析 11-2.2 讲过 Loadings()

  # (A2) 如果传入矩阵非空
  if (!IsMatrixEmpty(x = loadings)) {
  	#(B1) 基因数取最小值( 基因数，基因权重行数)
    nfeatures <- min(nfeatures, nrow(x = loadings))
    #(B2) 如果权重矩阵列数为0，则警告: 还没有投影，然后projected <- F，再载入 pca@feature.loadings
    if (ncol(x = loadings) == 0) {
      warning("Dimensions have not been projected. Setting projected = FALSE")
      projected <- FALSE
      loadings <- Loadings(object = x, projected = projected)
    }
    
    #(B3) 如果 dims 最小值 大于 权重矩阵列数，则报错
    if (min(dims) > ncol(x = loadings)) {
      stop("Cannot print dimensions greater than computed")
    }
    
    #(B4) 如果 dims 最大值 大于权重矩阵列数，则警告
    if (max(dims) > ncol(x = loadings)) {
      warning("Only ", ncol(x = loadings), " dimensions have been computed.")
      # 求交集: dims 和 每列编号
      # seq_len(length.out =3) #[1] 1 2 3
      dims <- intersect(x = dims, y = seq_len(length.out = ncol(x = loadings)))
    }

    #(B5) 遍历 dims
    for (dim in dims) {
      # features <- TopFeatures(
      #   object = x,
      #   dim = dim,
      #   nfeatures = nfeatures * 2,
      #   projected = projected,
      #   balanced = TRUE
      # )
      #(C1) 获取靠前的基因
      features <- Top(
        data = loadings[, dim, drop = FALSE], #取该列，保持df结构
        num = nfeatures * 2, #首尾分别取 n 个
        balanced = TRUE #正负 基因数相等，返回一个list
      )

      #"PC_1"
      cat(Key(object = x), dim, '\n')
      # 
      pos.features <- split(
        x = features$positive,
        f = ceiling(x = seq_along(along.with = features$positive) / 10)
      )

      cat("Positive: ", paste(pos.features[[1]], collapse = ", "), '\n')

      pos.features[[1]] <- NULL

      if (length(x = pos.features) > 0) {
        for (i in pos.features) {
          cat("\t  ", paste(i, collapse = ", "), '\n')
        }
      }

      neg.features <- split(
        x = features$negative,
        f = ceiling(x = seq_along(along.with = features$negative) / 10)
      )
      
      cat("Negative: ", paste(neg.features[[1]], collapse = ", "), '\n')
      neg.features[[1]] <- NULL
      if (length(x = neg.features) > 0) {
        for (i in neg.features) {
          cat("\t  ", paste(i, collapse = ", "), '\n')
        }
      }
    }
  }

  #(A3) 返回不可见对象
  return(invisible(x = x))
}














() UpdateKey()
./seurat-4.1.0/R/objects.R:2790:UpdateKey <- function(key) {
./seurat-object-4.0.4/R/utils.R:1018:UpdateKey <- function(key) {

#' Update a Key
#'
#' @param key A character to become a Seurat Key
#'
#' @return An updated Key that's valid for Seurat
#'
#' @section \code{Seurat} Object Keys:
#' blah
#'
#' @keywords internal
#'
#' @noRd
#'
UpdateKey <- function(key) {
  # 如果满足要求，直接返回: 字母数字 最后_结尾
  if (grepl(pattern = '^[[:alnum:]]+_$', x = key)) {
    return(key)
  # 如果不满足条件
  } else {
    # 则匹配出 字母数字 部分
    new.key <- regmatches(
      x = key,
      m = gregexpr(pattern = '[[:alnum:]]+', text = key)
    )
    # 用空字符 无缝拼接起来，最后添加下划线_
    new.key <- paste0(paste(unlist(x = new.key), collapse = ''), '_')
    # 如果只有下划线，则自动生成3个字符 后加下划线
    if (new.key == '_') {
      new.key <- paste0(RandomName(length = 3), '_')
    }
    # 警告: 1个或多个字符，并_结尾
    warning(
      "Keys should be one or more alphanumeric characters followed by an underscore, setting key from ",
      key,
      " to ",
      new.key,
      call. = FALSE,
      immediate. = TRUE
    )
    # 返回新字符串
    return(new.key)
  }
}














() DimHeatmap()
./seurat-4.1.0/R/visualization.R:43:DimHeatmap <- function(

#' Dimensional reduction heatmap
#'
#' Draws a heatmap focusing on a principal component. Both cells and genes are sorted by their
#' principal component scores. Allows for nice visualization of sources of heterogeneity in the dataset.
#' 画PC的热图。基因和细胞都按照PC打分排序。很好的可视化数据的变异来源。
#'
#' @inheritParams DoHeatmap
#' @param dims Dimensions to plot
#' @param nfeatures Number of genes to plot
#' @param cells A list of cells to plot. If numeric, just plots the top cells.
#' @param reduction Which dimensional reduction to use
#' @param balanced Plot an equal number of genes with both + and - scores.
#' @param projected Use the full projected dimensional reduction
#' @param ncol Number of columns to plot
#' @param fast If true, use \code{image} to generate plots; faster than using ggplot2, but not customizable
#' @param assays A vector of assays to pull data from
#' @param combine Combine plots into a single \code{\link[patchwork]{patchwork}ed}
#' ggplot object. If \code{FALSE}, return a list of ggplot objects
#'
#' @return No return value by default. If using fast = FALSE, will return a
#' \code{\link[patchwork]{patchwork}ed} ggplot object if combine = TRUE, otherwise
#' returns a list of ggplot objects
#'
#' @importFrom patchwork wrap_plots
#' @export
#' @concept visualization
#'
#' @seealso \code{\link[graphics]{image}} \code{\link[ggplot2]{geom_raster}}
#'
#' @examples
#' data("pbmc_small")
#' DimHeatmap(object = pbmc_small)
#'
DimHeatmap <- function(
  object,
  dims = 1, #作者的参数命名规则比较严谨：dims 表示可以是复数 1:5，而 ndims 就是一个数字。
  nfeatures = 30,
  cells = NULL,
  reduction = 'pca',
  disp.min = -2.5,
  disp.max = NULL,
  balanced = TRUE,
  projected = FALSE,
  ncol = NULL,
  fast = TRUE, #fast 为F才有返回值
  raster = TRUE, #栅格化
  slot = 'scale.data',
  assays = NULL,
  combine = TRUE
) {
  # (A1) 列数。如果为空，则看dims的长度，超过2的，按3列显示
  ncol <- ncol %||% ifelse(test = length(x = dims) > 2, yes = 3, no = length(x = dims))
  # (A2) 创建空向量，其实是一个list
  # > vector(mode = 'list', length = 2)
  #[[1]]
  #NULL
  #
  #[[2]]
  #NULL
  plots <- vector(mode = 'list', length = length(x = dims))

  #(A3) 默认 assays 是 "RNA"
  assays <- assays %||% DefaultAssay(object = object)

  #(A4) 设置数值的最大值。如果为空，则看slot，是'scale.data'的取2.5，否则(可能是data)取6
  disp.max <- disp.max %||% ifelse(
    test = slot == 'scale.data',
    yes = 2.5,
    no = 6
  )

  #(A5) 如果降维对象中记录的默认实验 不在 assays 中，则警告
  if (!DefaultAssay(object = object[[reduction]]) %in% assays) {
    warning("The original assay that the reduction was computed on is different than the assay specified")
  }



  ########### 开始排序 cell

  #(A6) 细胞，如果为空则使用全部细胞
  cells <- cells %||% ncol(x = object)

  #(A7) 如果 cells 是数字，则取每个pc的前n个 cid
  if (is.numeric(x = cells)) {
    cells <- lapply(
      X = dims, #对 dims 遍历
      FUN = function(x) {
        cells <- TopCells( #取每个 dim 的最高的细胞
          object = object[[reduction]], #传入降维类, 'pca'
          dim = x, #某个dim
          ncells = cells,
          balanced = balanced
        )
        # 如果正负平衡，则把负的细胞翻转顺序。
        # 原来是 c4,c3，翻转后是 c3, c4
        if (balanced) {
          cells$negative <- rev(x = cells$negative)
        }
        # 然后解开list：类似 c1-c4这样的字符串
        cells <- unlist(x = unname(obj = cells))
        return(cells)
      }
    ) #lapply返回值是list
  }
  #> lapply(X=1:3, FUN=function(x){
  #+     return(paste0("Cell",1:x))
  #+ })
  #[[1]]
  #[1] "Cell1"
  #
  #[[2]]
  #[1] "Cell1" "Cell2"
  #
  #[[3]]
  #[1] "Cell1" "Cell2" "Cell3"

  # (A8) 如果 cells 不是 list形式，则生成list，键值对数量和dims相等
  if (!is.list(x = cells)) {
    cells <- lapply(X = 1:length(x = dims), FUN = function(x) {return(cells)})
  }



  ########### 开始排序 gene

  #(A9) 找每个 dim 的load最高的基因
  features <- lapply(
    X = dims,
    FUN = TopFeatures,
    # 后面几个是 TopFeatures() 的其他几个参数
    object = object[[reduction]],
    nfeatures = nfeatures,
    balanced = balanced,
    projected = projected
  )

  # (A10) 解开list，取唯一值，获得全部基因名
  features.all <- unique(x = unlist(x = features))

  # (A11) 如果 assays 长度大于1
  if (length(x = assays) > 1) {
    # 
    features.keyed <- lapply(
      X = assays, # 对 assays 遍历
      FUN = function(assay) {
        # 求基因的交集
        features <- features.all[features.all %in% rownames(x = object[[assay]])]
        # 如果长度大于0，则返回 key 前缀 + 基因名
        if (length(x = features) > 0) {
          return(paste0(Key(object = object[[assay]]), features))
        }
      }
    )
    # 过滤掉空值。为什么会有空值? 
    features.keyed <- Filter(f = Negate(f = is.null), x = features.keyed)
    # 解开list，获取带有 key前缀的基因名列表
    features.keyed <- unlist(x = features.keyed)
  # 如果 assays 长度为 1
  } else {
    # 免去前缀
    features.keyed <- features.all
    # 默认实验为 assays，是不是循环赋值了？
    DefaultAssay(object = object) <- assays
  }

  #(A12) 获取数据, 该函数在 源码解析9-2.1 讲过
  data.all <- FetchData(
    object = object, #Seurat 对象
    vars = features.keyed, #基因名字，对多个 assays的，带有 其key前缀。
    cells = unique(x = unlist(x = cells)), #解开list，获取细胞id
    slot = slot #默认是 'scale.data'
  )

  #(A13) 按照设置的最值截断数据
  data.all <- MinMax(data = data.all, min = disp.min, max = disp.max)
  # (A14) 数据的范围
  data.limits <- c(min(data.all), max(data.all))


  # if (check.plot && any(c(length(x = features.keyed), length(x = cells[[1]])) > 700)) {
  #   choice <- menu(c("Continue with plotting", "Quit"), title = "Plot(s) requested will likely take a while to plot.")
  #   if (choice != 1) {
  #     return(invisible(x = NULL))
  #   }
  # }
  # 如果快速模式，则重新定义行列数
  if (fast) {
    nrow <- floor(x = length(x = dims) / 3.01) + 1
    orig.par <- par()$mfrow
    par(mfrow = c(nrow, ncol))
  }

  #(A15) 遍历 dim
  for (i in 1:length(x = dims)) {
    # 取出每个 features list中的基因，两层list，第一层是dim，第二层是 pos / neg
    # 相当于 x6 x5,  x2, x1; 原来是降序排列，现在倒过来，是升序排列了
    dim.features <- c(features[[i]][[2]], rev(x = features[[i]][[1]]))
    # (B2)匹配出 带有key 前缀的基因名字
    dim.features <- rev(x = unlist(x = lapply(
      X = dim.features, #对每个基因 遍历
      FUN = function(feat) { # 匹配 (带key前缀的)features.keyed 中结尾是该基因的条目，返回匹配值。
        return(grep(pattern = paste0(feat, '$'), x = features.keyed, value = TRUE))
      }
    )))
    # (B3)属于该 dim 的cell id
    dim.cells <- cells[[i]]
    # (B4) 从表达值df中取子集
    data.plot <- data.all[dim.cells, dim.features]
    
    #(B5) 如果fast==T，则使用 SingleImageMap() 函数
    if (fast) {
      SingleImageMap( #这个竟然没有返回值？竟然还是默认值！！
        data = data.plot,
        title = paste0(Key(object = object[[reduction]]), dims[i]), #"PC_1"
        order = dim.cells
      )
    # 如果 fast==F，则返回的是 SingleRasterMap() 函数
    } else {
      plots[[i]] <- SingleRasterMap(
        data = data.plot,
        raster = raster, #默认是 栅格化
        limits = data.limits, #数据的2个最值
        cell.order = dim.cells,
        feature.order = dim.features
      )
    }
  }


  #(A16) 如果 fast==T，则复原图形参数，并返回不可见的NULL
  if (fast) {
    par(mfrow = orig.par)
    return(invisible(x = NULL)) #默认这里返回了
  }

  # (A17) 如果要拼接，则拼接并返回
  if (combine) {
    plots <- wrap_plots(plots, ncol = ncol, guides = "collect")
  }
  return(plots)
}










() Top()
./seurat-object-4.0.4/R/utils.R:912:Top <- function(data, num = 20, balanced = FALSE) {
./seurat-4.1.0/R/objects.R:2653:Top <- function(data, num, balanced) {


#' Get the top
#'
#' @param data Data to pull the top from
#' @param num Pull top \code{num}
#' @param balanced Pull even amounts of from positive and negative values
#'
#' @return The top \code{num}
#'
#' @importFrom utils head tail
#'
#' @keywords internal
#'
#' @noRd
#'
Top <- function(data, num = 20, balanced = FALSE) {
  # 行数
  nr <- nrow(x = data)
  # 如果要打印的基因数超过总行数，则警告，并设置打印的行数为 总行数
  if (num > nr) {
    warning(
      "Requested number is larger than the number of available items (",
      nr,
      "). Setting to ",
      nr ,
      ".",
      call. = FALSE
    )
    num <- nr
  }
  # 如果总行数为1，则不用平衡基因数了
  balanced <- ifelse(test = nr == 1, yes = FALSE, no = balanced)

  # 按1列矩阵，选取前n个基因
  top <- if (isTRUE(x = balanced)) { #默认走这里
    # 分正负，所以除以2
    num <- round(x = num / 2)
    # 降序排列，最大的在前面
    data <- data[order(data, decreasing = TRUE), , drop = FALSE]

    # 取前n个
    positive <- head(x = rownames(x = data), n = num)
    # 取后n个，并倒序
    negative <- rev(x = tail(x = rownames(x = data), n = num))

    # 作者想在这里去重复，不过这样写肯定不对。 //bug: 记录到本文，不提交PR了。
    # 因为99.99%的情况用不到，单细胞的基因数至少上千个，而每个PC要打印出来的不超过100个，正负top100都没有交集。
    # remove duplicates
    if (positive[num] == negative[num]) {
      negative <- negative[-num]
    }
    # 返回一个list
    list(positive = positive, negative = negative)
  } else { #默认是不走这里
    # 按绝对值降序排列
    data <- data[rev(x = order(abs(x = data))), , drop = FALSE]
    # 选前n个
    top <- head(x = rownames(x = data), n = num)
    # 在按算术值排序
    top[order(data[top, ])]
  }
  return(top)
}










() TopCells()
./seurat-4.1.0/R/objects.R:753:TopCells <- function(object, dim = 1, ncells = 20, balanced = FALSE, ...) {

#' Find cells with highest scores for a given dimensional reduction technique
#'
#' Return a list of genes with the strongest contribution to a set of components
#'
#' @param object DimReduc object
#' @param dim Dimension to use
#' @param ncells Number of cells to return
#' @param balanced Return an equal number of cells with both + and - scores.
#' @param ... Extra parameters passed to \code{\link{Embeddings}}
#'
#' @return Returns a vector of cells
#'
#' @export
#' @concept objects
#'
#' @examples
#' data("pbmc_small")
#' pbmc_small
#' head(TopCells(object = pbmc_small[["pca"]]))
#' # Can specify which dimension and how many cells to return
#' TopCells(object = pbmc_small[["pca"]], dim = 2, ncells = 5)
#'
TopCells <- function(object, dim = 1, ncells = 20, balanced = FALSE, ...) {
  # dim(Embeddings(pbmc@reductions$pca))
  # Embeddings(pbmc@reductions$pca)[1:2,1:3]
  #                        PC_1       PC_2       PC_3
  # AAACATACAACCAC-1 -4.7298963 -0.5182652 -0.7809100
  # AAACATTGAGCTAC-1 -0.5176254  4.5923068  5.9605692
  embeddings <- Embeddings(object = object, ...)[, dim, drop = FALSE] #见 源码解析 3-1.3
  return(Top(
    data = embeddings,
    num = ncells,
    balanced = balanced
  ))
}




() TopFeatures()
./seurat-4.1.0/R/objects.R:715:TopFeatures <- function(

#' Find features with highest scores for a given dimensional reduction technique
#'
#' Return a list of features with the strongest contribution to a set of components
#'
#' @param object DimReduc object
#' @param dim Dimension to use
#' @param nfeatures Number of features to return
#' @param projected Use the projected feature loadings
#' @param balanced Return an equal number of features with both + and - scores.
#' @param ... Extra parameters passed to \code{\link{Loadings}}
#'
#' @return Returns a vector of features
#'
#' @export
#' @concept objects
#'
#' @examples
#' data("pbmc_small")
#' pbmc_small
#' TopFeatures(object = pbmc_small[["pca"]], dim = 1)
#' # After projection:
#' TopFeatures(object = pbmc_small[["pca"]], dim = 1,  projected = TRUE)
#'
TopFeatures <- function(
  object,
  dim = 1,
  nfeatures = 20,
  projected = FALSE,
  balanced = FALSE,
  ...
) {
  # 见 源码解析 11-2.2
  loadings <- Loadings(object = object, projected = projected, ...)[, dim, drop = FALSE]
  return(Top(
    data = loadings,
    num = nfeatures,
    balanced = balanced
  ))
}












() MinMax()
./seurat-4.1.0/R/utilities.R:1100:MinMax <- function(data, min, max) {

#' Apply a ceiling and floor to all values in a matrix
#'
#' @param data Matrix or data frame
#' @param min all values below this min value will be replaced with min
#' @param max all values above this max value will be replaced with max
#' @return Returns matrix after performing these floor and ceil operations
#' @export
#' @concept utilities
#'
#' @examples
#' mat <- matrix(data = rbinom(n = 25, size = 20, prob = 0.2 ), nrow = 5)
#' mat
#' MinMax(data = mat, min = 4, max = 5)
#'
MinMax <- function(data, min, max) {
  data2 <- data
  data2[data2 > max] <- max
  data2[data2 < min] <- min
  return(data2)
}









() SingleImageMap () 默认
./seurat-4.1.0/R/visualization.R:7412:SingleImageMap <- function(data, order = NULL, title = NULL) {

使用 base-R 的画图命令 image() 画热图。

#' A single heatmap from base R using \code{\link[graphics]{image}}
#'
#' @param data matrix of data to plot
#' @param order optional vector of cell names to specify order in plot
#' @param title Title for plot
#'
#' @return No return, generates a base-R heatmap using \code{\link[graphics]{image}}
#'
#' @importFrom graphics axis image par plot.new title
#'
#' @keywords internal
#'
#' @export
#'
SingleImageMap <- function(data, order = NULL, title = NULL) {
  # 如果order 非空，则按行对data排序
  if (!is.null(x = order)) {
    data <- data[order, ]
  }
  # 单位line 设置边距，底部起，顺时针 c(bottom, left, top, right)
  par(mar = c(1, 1, 3, 3))
  # 新图形
  plot.new()
  # image() 是画热图的，输入矩阵
  image(
    x = as.matrix(x = data), #必须输入矩阵，df不行
    axes = FALSE,
    add = TRUE, #是否加到现有图形上
    col = PurpleAndYellow() #颜色，应该是一个渐变色
  )

  # 添加坐标轴
  axis(
    side = 4, #1=below, 2=left, 3=above and 4=right.
    at = seq(from = 0, to = 1, length = ncol(x = data)), #位置是 0-1 之间，几列分成几份
    labels = colnames(x = data), #刻度的标签是 列名
    las = 1, #总是水平
    tick = FALSE, #不显示刻度线
    mgp = c(0, -0.5, 0), #3个距离
    cex.axis = 0.75
  )
  # 图片的标题
  title(main = title)
}










() SingleRasterMap
./seurat-4.1.0/R/visualization.R:7479:SingleRasterMap <- function(
使用 ggplot2 画热图，默认栅格化。


#' A single heatmap from ggplot2 using geom_raster
#'
#' @param data A matrix or data frame with data to plot
#' @param raster switch between geom_raster and geom_tile #切换2个画图函数
#' @param cell.order ...
#' @param feature.order ...
#' @param colors A vector of colors to use
#' @param disp.min Minimum display value (all values below are clipped)
#' @param disp.max Maximum display value (all values above are clipped)
#' @param limits A two-length numeric vector with the limits for colors on the plot
#' @param group.by A vector to group cells by, should be one grouping identity per cell
#'
#' @return A ggplot2 object
#
#' @importFrom ggplot2 ggplot aes_string geom_raster scale_fill_gradient
#' scale_fill_gradientn theme element_blank labs geom_point guides
#' guide_legend geom_tile
#'
#' @keywords internal
#'
#' @export
#
SingleRasterMap <- function(
  data,
  raster = TRUE,
  cell.order = NULL,
  feature.order = NULL,
  colors = PurpleAndYellow(),
  disp.min = -2.5,
  disp.max = 2.5,
  limits = NULL,
  group.by = NULL
) {
  # (A1) 按上下极值截断数据
  data <- MinMax(data = data, min = disp.min, max = disp.max)
  
  # (A2) 宽变长，共3列: rows         cols vals
  data <- Melt(x = t(x = data))
  # 重命名 列名
  colnames(x = data) <- c('Feature', 'Cell', 'Expression')
  
  # (A3) 整理数据列
  # 如果基因名排序参数非空，则基因列 强转为 因子
  if (!is.null(x = feature.order)) {
    data$Feature <- factor(x = data$Feature, levels = unique(x = feature.order))
  }

  # 如果 细胞名排序参数非空，则细胞列 强转为 因子
  if (!is.null(x = cell.order)) {
    data$Cell <- factor(x = data$Cell, levels = unique(x = cell.order))
  }

  # 如果 分组 参数非空，则添加一列
  if (!is.null(x = group.by)) {
    data$Identity <- group.by[data$Cell] #这个group.by 是不是应该是 named vector
    # 例:
    #> a1=1:10; names(a1)=c( paste0("X",seq(10,5,-1)), paste0("X",seq(1,4,1)) )
    #> a1
    # X10  X9  X8  X7  X6  X5  X1  X2  X3  X4 
    #   1   2   3   4   5   6   7   8   9  10 
    #> a1[ paste0("X", 1:10)]
    # X1  X2  X3  X4  X5  X6  X7  X8  X9 X10 
    #  7   8   9  10   6   5   4   3   2   1 
  }

  # (A4) 定义表达数据的上下限
  limits <- limits %||% c(min(data$Expression), max(data$Expression))
  # 如果界限长度不是2，或者不是数字，则报错。
  if (length(x = limits) != 2 || !is.numeric(x = limits)) {
    stop("limits' must be a two-length numeric vector")
  }


  ##################
  # 开始画图
  ##################
  # (A5) 画热图的函数，默认T，选择 geom_raster
  my_geom <- ifelse(test = raster, yes = geom_raster, no = geom_tile)

  # 画图
  plot <- ggplot(data = data) +
    my_geom(mapping = aes_string(x = 'Cell', y = 'Feature', fill = 'Expression')) +

    theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) + #去掉x轴文字、刻度
    scale_fill_gradientn(limits = limits, colors = colors, na.value = "white") + #数值范围，渐变色
    labs(x = NULL, y = NULL, fill = group.by %iff% 'Expression') + #标签
    WhiteBackground() + NoAxes(keep.text = TRUE) #白色主题，没有坐标轴


  # (A6) 如果  group.by 非空
  if (!is.null(x = group.by)) {
    # 散点图，全透明，就是为了获得一个图例！?
    plot <- plot + geom_point(
      mapping = aes_string(x = 'Cell', y = 'Feature', color = 'Identity'),
      alpha = 0
    ) +
      guides(color = guide_legend(override.aes = list(alpha = 1))) #一个 color的图例，代替fill的图例
  }

  return(plot)
}










() Melt()
./seurat-4.1.0/R/utilities.R:2122:Melt <- function(x) {

见本文 3.11





() WhiteBackground()

见本文 3.12


() NoAxes()
见本文 3.12

